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One of the main features of astrophysical shocks is their ability to acceler- 
ate particles to extremely high energies. The leading acceleration mechanism, 
the diffusive shock acceleration is reviewed. It is demonstrated that its effi- 
ciency critically depends on the injection of thermal plasma into acceleration 
which takes place at the subshock of the collisionless shock structure that, in 
turn, can be significantly smoothed by energetic particles. Furthermore, their 
inhomogeneous distribution provides free energy for MHD turbulence regulat- 
ing the subshock strength and injection rate. Moreover, the MHD turbulence 
confines particles to the shock front controlling their maximum energy and 
bootstrapping acceleration. Therefore, the study of the MHD turbulence in a 
compressive plasma flow near a shock is a key to understanding of the entire 
process. The calculation of the injection rate became part of the collisionless 
shock theory. It is argued that the further progress in diffusive shock acceler- 
ation theory is impossible without a significant advance in these two areas of 
plasma physics. 
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I. INTRODUCTION 



Over the last few years new observations and missions, e.g., Energetic Gamma-ray Ex- 
periment Telescope (EGRET) Chandra p[, TeV- astronomy 0, revolutionized the mea- 
surement of radiation from a variety of objects in the Universe. In most cases the primary 
source of the radiation is believed to be accelerated charged particles, often of remarkably 
high energies, such as lO^^eV or even higher [Q, usually referred to as the ultra high en- 
ergy cosmic rays (UHECR). The accelerated particles themselves (we will also use the term 
cosmic rays, CRs) are in many cases generated by shock waves (shocks). Note, that the 
latter are the major events where the huge energy of stars, Supernovae's (SN) or Black 
holes is released in bulk gas motions and ultimately dissipated. The most successful particle 
acceleration mechanism is, perhaps the diffusive shock acceleration (DSA) which is a 
variant of the original Fermi idea (1949), also known as the First order Fermi acceleration 
process. According to this mechanism particles gain energy by bouncing off hydromagnetic 
disturbances frozen into the converging upstream and downstream flow regions near a shock. 
Clearly, the understanding of this mechanism is critical for radiation models since the pri- 
mary particle spectrum is one of their most important "input" characteristics. The success 
of this mechanism has been due to the following appealing features 

1. it reproduces the power-law energy spectra with an index remarkably similar to that 
inferred from the observations of the galactic cosmic rays (CRs) [^,^. 

2. convincing, direct observational evidence of its operation at interplanetary shocks and 
the earth bow shock 1^-0 

3. the absence of any obvious intrinsic limitation to the maximum particle energy achiev- 
able by this process 

Often, however, this mechanism is applied in its simplified version, the so called test par- 
ticle (TP) or linear approximation (that neglects the backreaction of accelerated particles 
on the shock structure and produces a simple E~'^ particle energy distribution). Although 
more realistic, some nonlinear theories do exist, they suffer from parameterization of par- 
ticle transport coefficients such as the pitch-angle scattering, spatial diffusivity and plasma 
heating. 

Also, it should be emphasized that the statement (2) above is only a prima facie evi- 
dence for the responsibility of this mechanism for galactic CRs, alluded to in (1). Despite 
similarity between physical parameters 0, the life-times and extensions of shocks in these 
two environments are vastly different (up to a factor ~ 10^''). This requires different ap- 
proximations for these shocks. Namely, the treatment of large astrophysical shocks such as 
Supernova remnant (SNR) shells, shocks in the lobes of radio galaxies or even larger shocks 
in clusters of galaxies must include the back reaction of accelerated particles on the shock 
structure, and thus be intrinsically nonlinear. This, in turn, necessitates the selfconsistent 
treatment of the above mentioned anomalous transport phenomena. 

From the perspective of the widely accepted SNR shock hypothesis of the origin of 
galactic CRs, perhaps, the most critical impact of this treatment would be on the statement 
(3) above. It implies that the maximum energy is limited only by the time available for 
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acceleration and by the size of accelerating object (roughly not to be exceeded by the diffusive 
escape length for particles, which grows with energy). This is becoming a "hot" issue 
in view of the lack of evidence for TeV-protons in SNR shocks (despite the signature of 
GeV-protons, EGRET), which might indicate that the particle spectrum either cuts off 
somewhere between GeV and TeV energies (i.e., in a currently uncovered by any instrument 
energy range) or probably significantly steepens there. We must await observations from 
Gamma-ray Large Area Space Telescope (GLAST) and the next generation of the imaging 
atmospheric Cherenkov telescopes that are currently being built to see where and how 
the spectrum disappears or changes. To understand how this may happen, one needs to 
consider all the crucial requirements for this mechanism, which are 



• the ability of some fraction of thermal particles downstream to make the first step to 
acceleration, i.e., to return upstream (injection problem) 

• good confinement of accelerated particles near the shock front (to continue energy 
gain) 

• the ability of the shock to withstand the pressure of accelerated particles (the shock 
smoothing, or shock robustness problem) 

These requirements are also associated with the main difficulties of the theory and are, in 
fact, strongly related. We consider them briefly, in order. 



A. Injection problem 



Injection is the process of initial particle energization. Within the widely accepted "ther- 
mal leakage" scenario of injection, the way protons enter the shock acceleration process is 
physically the same as when they are accelerated afterwards (see, e.g., Ref. [0). After 
thermalization downstream, certain fraction of particles will catch up with the shock. Their 
further leakage upstream generates Alfven waves via the cyclotron resonance kpfi ^ eB/c 
where /i is the cosine of the proton pitch angle, and p is the particle momentum. These waves 
do two things. First, they self-regulate the thermal particle leakage (injection) by trapping 
particles downstream when their leakage becomes too strong and therefore the wave ampli- 



tude too large as well [|TT|. Second, the large amplitude waves scatter already accelerating 
particles in pitch angle, thus ensuring their diffusive confinement near the shock. 

According to the above resonance condition, waves excited by protons (with kpp ~ 1) 
are too long to scatter electrons, so that they interact with them adiabatically {kpe -C 
1). Therefore, electrons need a separate injection scenario. One suggestion is that they 
might be injected via scattering on self-generated whistler waves [0. According to another 



mechanism suggested in Ref. []T3[, electrons are extracted directly from the thermal pool 
near the shock front via their interaction with proton-generated lower-hybrid waves and 
electrostatic barrier formed there due to this interaction. 

It should be noted that when the both species become relativistic and, if the synchrotron 
cooling time is longer than the acceleration time, the electron spectrum should be identical 
to that of the protons. However, it is electrons for which we have convincing direct evidence 
of Fermi acceleration in SNR shells. There are well documented measurements of nonthermal 
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electron emission in radio, x rays and possibly also in 7 rays, presented e.g., in Refs. [jT4[, 
15| and |16|, respectively. Although electrons play perhaps dynamically unimportant role 
in the shock structure, they should trace the dynamically important proton spectra. 



B. Particle Confinement near the shock 

As in the injection phase, further particle confinement near the shock is supported by 
self-generated Alfven wave turbulence, since the energetic particle distribution ahead of the 
shock is ion-cyclotron unstable (see, e.g., Ref. pj). Due to rapid pitch-angle scattering on 
these waves, particles cross the shock repeatedly, thus gaining energy. Clearly, the scattering 
frequency (turbulence level) sets up the acceleration rate. Note that the existing background 
magneto-hydrodjTiamic (MHD) turbulence in the interstellar medium (ISM) would support 
only very slow acceleration. A long-standing problem is that the wave generation process 
is likely to be so robust that wave amplitude may far exceed the level admissible by the 
quasi-linear theory. It is usually expected that the turbulence saturates at a level 6B ~ Bq. 
Thus, particle scattering must occur via strongly nonlinear wave-wave and wave-particle 
interactions so that the conventional quasi-linear description should be replaced by a non- 
perturbative approach. 

MHD simulations significantly advanced recently (see, eg |]17[ and references therein) but 



they would not suffice alone to self-consistently describe the wave generation by turbulently 
confined particles. The main difficulty is, as it will be seen in the sequel, the enormous ex- 
tension of particle and wave spectra. It seems desirable to combine simulations in restricted 
but critical parts of the phase space (as the short-wave, low-energy part that controUs the 
injection and the long-wave, high energy part where particle losses occur) with an analytical 
approach in extended but tractable parts of the phase space where particle and wave spectra 



and even the fiow profile exhibit relatively simple, scale invariant behaviour . 

It should be evident that a purely numerical approach would encounter serious difficulties 
if applied to the shock acceleration of UHECRs. Indeed, even if we had numerical solution 
of the problem of CR acceleration in SNR shocks, i.e., up to energies ~ 10^^ eV, we would 
need to extend its dynamical range by a factor of ~ 10^. Since the integration domain in 
configuartion space is typically proportional to the maximum energy, this would mean ~ 10^° 
times larger phase space. It seems natural to use this large size of the phase space instead 
of letting computer to process decade after decade of it with basically the same physics. 
The analytic approach that we introduce in Sec. [il B| below utilizes the corresponding small 



parameters and therefore appears particularly suitable to the problem of acceleration of 
UHECRs. 



C. Shock robustness 

Due to the secular growth of the pressure of accelerated particles, the linear solution 
becomes invalid for shocks of sufficient life time, scale L and magnetic field B (the most 
critical parameter is BL since the maximum particle energy scales as Emax ~ {u/c)eBL, 
provided that enough time is available for acceleration, as discussed, e.g., in Ref. 0). The 
main nonlinear effect is due to the back reaction of energetic particles on the fiow which 
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reduces the velocity jump at the flow discontinuity (subshock) through the deceleration 
and heating of the plasma in front of it (i.e., in the so called CR precursor), Fig.|l[ The 
total shock compression, however, increases due to the decrease in the adiabatic index 7 
caused by the presence of relativistic particles and, even more importantly, by their escape 
from the system (akin to radiative shocks). At the same time the maintenance of a finite 
subshock within the global shock structure (subshock itself, plus CR precursor) is critical 
to the injection process and thus for acceleration in general. A mathematical limitation 
of the TP approximation can be obtained in terms of particle energy and injection rate 
u, namely, uy/Emax ~ 1 which states that the pressure of accelerated particles must 



remain smaller than the shock ram pressure. It is also clear from this condition that the 
above three issues are strongly coupled. Injection (z/) and particle confinement (which sets 
Emax) determine the shock structure and are, in turn, regulated by it. The acceleration 
time scale Tacc (or particle diffusivity n) and the precursor turbulent heating rate are also 
implicitly involved in this "feedback loop". There are indeed too many shock variables to 
obtain reliable prediction by simply scanning the parameter space. Therefore, they need to 
be calculated self-consistently before the particle spectra can be calculated and compared 
with the observations. 



II. NONLINEAR THEORIES 



A typical nonlinearly accelerating shock is illustrated in Fig.|l]. The most "visible" non- 
linear effects are: (1) the deceleration of the flow upstream (see the flow profile on the top 
of Fig. |l|) by the pressure of accelerated particles, (2) subshock reduction, as a consequence 
(its strength uo/u2 may become small compared to the total compression Ui/u2 (3) bending 
of the B-field due to the compression of its tangential component because of the frozen in 
condition, uBt = const, while the normal component is conserved due to divB = 0. 



A. Two-fluid model 



Initially, the back reaction of accelerated particles (CRs) on the shock structure was 
studied within the two- fluid model (TFM). This model treats CRs as a second fluid carrying 
the momentum and energy across the shock, but not the mass. Complete solutions given in 
Refs. [pOtpT| indeed revealed a very strong back reaction of the CRs onto the bulk plasma 



flow, leading to a bifurcation of the simple linear solution into a strongly nonlinear (efficient) 
solution with the acceleration efficiency approaching (in strong shocks) 100% due to the 
nonlinearly enhanced shock compression. This gives rise to a formally diverging (i.e., in 
reality strongly dependent on the maximum energy) CR pressure which makes the particle 
losses at the highest energies as important dynamically as injection, and strongly related to 
it through their feedback on the subshock strength. Moreover, these losses are controlled by 
the Alfven turbulence which, in turn, depends on the CR distribution in the shock precursor 
and thus on the losses themselves. Unfortunately, the TFM, being a hydrodynamic theory, 
cannot be closed in such a way that these essentially kinetic effects are properly represented. 
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B. Kinetic Theory 



The earlier kinetic theories also demonstrated that upon accumulating enough CR energy, 
strong shocks develop into the nonlinear regime The minimal kinetic theory that captures 
bifurcations of the DSA can be formulated as follows. 

We describe the distribution of high energy particles (CRs) by the so called diffusion- 
convection equation (see, e.g., Refs. f^J^). The gaseous discontinuity (the subshock) is 
assumed to be located at x = and, for convenience, we flip the x- coordinate in Fig.^, so 
that the upstream side is x > half-space. Thus, the flow velocity in the shock frame can be 
represented as V{x) = ~u{x) where the (positive) flow speed u{x) jumps from U2 = w(0— ) 
downstream to uq = u{0+) > U2 across the subshock and then gradually grows upstream 
up to Ml = ^(+00) > Uq. In a steady state the equation can be written as 

+ nip)^ = -—p^ (1) 
dx dx'^ 3dx dp'' 

where f{x,p) is the isotropic (in the local fluid frame) part of the particle distribution. 
This is assumed to vanish far upstream {f ^ 0, x —>■ 00), while the only bounded solution 
downstream is obviously f{x,p) = fo{p) = f{0,p)- The most plausible assumption about the 
cosmic ray diffusivity k{p) is that of the Bohm type, i.e., n{p) = Kp^ / ^J1 + (the particle 
momentum p is normalized to mc). In other words k scales as the gyroradius, k ~ Tg{p)- The 
reference diffusivity K depends on the 6B/B level of the Alfvenic turbulence that scatters 
the particles in pitch angle. The minimum value for K would be ~ mc^/cB if 6B ~ B. 
Note that the replacement of this plain parameterization by a selfconsistent solution for the 
spectrum of turbulence driven by an inhomogeneous distribution of accelerated particles is a 
challenge to plasma physics. The existing quasi-linear approaches clearly fail in the efficient 
acceleration regime due to unacceptable wave amplitudes. 

The determination of u{x) in eq.([l|) requires three further equations. The first one is the 
conservation of the momentum flux in the smooth part of the shock transition (x > 0, i.e., 
in the CR-precursor) 

+ pu"^ = p^uj, x>0 (2) 
where Pc is the pressure of the high energy particles 

Pc{x) = —mc / -==f{p,x) (3) 
Jpo VP + 1 

It is assumed here that there are no particles with momenta p > pi (they leave the shock 
vicinity because there are no MHD waves with sufficiently long wave length A, since the 
cyclotron resonance requires p ~ A). The momentum region < p < Pq cannot be described 
by equation ([1|) and the behavior of f{p) at p ~ po is described by the injection parameters 
Po and f{po) HHI]. The plasma density p{x) can be eliminated from equation (H) by using the 
continuity equation pu = piUi. Finally, the subshock strength r^ can be expressed through 
the Mach number M sX x = 00 

r = I±i (4) 

' U2 7 - 1 + 2i?^+iM-2 ^ ^ 
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where the precursor compression R = ui/uq and 7 is the adiabatic index of the plasma. 

The system of equations (|l|,|^J^) describes in a self-consistent manner the particle spec- 
trum and the flow structure, although under parameterization of such critical quantities as 
u and pi. Our poor knowledge of the maximum momentum pi is related to the prescribed 
form of particle diffusivity k{p). 

It is useful to reduce this system to one integral equation [|l^]. A key dependent variable 
is an integral transform of the flow profile u{x) with a kernel suggested by an asymptotic 
solution of the system (|l|)-(@) which has the form 



f{x,p) = /o(p)exp 



3k 



where 



u{x')dx' 



is the flow potential and the spectral index downstream q{p) = — (iln/o/cilnp. The integral 
transform is as follows 



1 f°° 
U{p) = — / exp 
Ui Jo- 



qjp) 

3k{p) 



(5) 



and it is related to q{p) through the following formula 

dlnU 3 



(6) 



The physical meaning of the function U{p) is very simple. It reflects the degree of shock 
modification. Namely, a function U (p) + U2 is an effective flow velocity upstream as seen 
by a particle with momentum p that diffusively escapes ahead of the shock to a point x 
where the flow speed is u{x) = U{p) +U2- Once U{p) is found, both the flow profile and the 
particle distribution can be determined by inverting transform (|^) and integrating equation 
(I). Now, using the linearity of equation (||) {pu = const), we derive the integral equation 
for U by applying the transformation @ to the x— derivative of equation [|l^. The 
result reads 



Uit) 



1 



+ 



Kpo J to 



dt' 



+ 



K{t') K{t)q{t) 



-1 



3 


fi' dt" ' 


Rvs J 


to U{t")_ 



where t = Inp, to,i = lnpo,i- Here the injection parameter 

47r mc^ 



(7) 



3 piu 



rPo/o(Po 



(8) 



is related to R by means of the following equation 



7 




to U[t) 



dt' 



Rr, Jto U{t') 



(9) 



The equations (^J^,^) form a closed system that can be easily solved, 
Asymptotic analytic solutions are also available 
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e.g., numerically. 
The physical quantities involved are: 
the far upstream Mach number M (which is given, external parameter); internal parameters: 
Mq (the Mach number at the subshock), the injection rate z/, the particle maximum momen- 
tum Pmax = Pi, and particle diffusivity n{p,x). These internal parameters must be deter- 
mined self-consistently, but currently are parameterized or calculated using some simplifying 
assumption. For example, assuming that there is no turbulent heating in the shock precursor 
(which is doubtful in such a turbulent environment), the parameter Mq = M/R^'^^^^^'^ [see 
eq.(l)]. 



III. CRITICAL NATURE OF ACCELERATION PROCESS 

The presence of bifurcation in this acceleration process is best seen in variables R, v. The 
quantity — 1 is a measure of shock modification produced by CRs, in fact {R — 1)/R = 
Pc{0)/piul [Eq.(^] and may be regarded as an order parameter. The injection rate u 
characterizes the CR density at the shock front and can be tentatively treated as a control 
parameter. It is convenient to plot the function z/(_R) instead of R{i^) [using equations (^ 
and i)], Fig.|. 

In fact, the injection rate u at the subshock should be calculated given rs{R) with the 
self-consistent determination of the flow compression R on the basis of the -R(z^) dependence 
obtained. However, in view of a very strong, even nonunique dependence -R(z/), this solution 
can be physically meaningful only in regimes far from criticality, i.e., when R ^ 1 (test par- 
ticle regime) or i? ^ 1 (efficient acceleration). There are, however, self-regulating processes 
that should drive the system towards critical region (where -R(z^) dependence is very sharp) 
18| , |22[| . First, if u is subcritical it will inevitably become supercritical when pi grows in 



course of acceleration. Once it happened, however, the strong subshock reduction (equation 
Q) will reduce u as well and drive the system back to the critical regime, Fig.|^. 

The maximum momentum pi is subject to self- regulation as well. Indeed, when i? 3> 1, 
the generation and propagation of Alfven waves is characterized by strong inclination of the 
characteristics of wave transport equation towards larger wavenumbers k on the k — x plane 
due to wave compression. Thus, considering particles with p ^ pi inside the precursor, one 
sees that they are in resonance with waves that must have been excited by particles with 
p > Pi further upstream but, there are no particles with p > pi. Therefore, the required 
waves can be excited only locally by the same particles with p ^ pi which substantially 
diminishes the amplitude of waves that are in resonance with particles from the interval 
Pi/ R < p < Pi. (The left inequality arises from the resonance condition kcp ~ eB/mc and 
the frequency conservation along the characteristics ku{x) = const). This will worsen the 
confinement of these particles to the shock front. The quantitative study of this process 
is another challenge to the theory of plasma turbulence. What can be inferred from Fig.|^ 
now, is that the decrease of pi straightens out and rise the curve z/(i?), so that it returns to 
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the monotonic behaviour. However once the actual injection becomes subcritical (and thus 
i? — 1) then pi will grow again restoring the two extrema on the curve z/(i?). 

Also the turbulent precursor heating straightens up the bifurcation diagram and returns it 
to the critical state. We illustrate this in Fig.^ for different heating efficiencies a, introduced 
phenomenologically through the decrease in the flow Mach number across the CR precursor 

= M-'R'/' + - l)Pmacc (10) 

Here the first term is due to the familiar adiabatic compression while the second describes the 
turbulent heating, assumed to be proportional to the CR pressure contrast oc i? — 1 (source 
of the free energy for turbulence) as well as to the precursor length oc Pmax- The plots in 
Figs. demonstrate that parameterization approaches in which the internal parameters 
(like a and Pmax in many simulations) are treated as given may be useful only when (and if) 
these parameters are well outside the critical region which is marked by very sharp or even 
nonunique parameter dependence. For example, as seen from Fig.^, if a (which is essentially 
unknown) lies within the interval 10^® — 10^^, then there are no means to calculate the flow 
structure (R), and thus the particle spectrum with a reasonable accuracy. To treat the 
injection rate in the critical region as a "control" parameter is equally useless. Figs. 0,^. 

It appears to be more productive to assume a marginal state in which the maximum and 
the minimum of z^(-R) merge, at least in a sense of an averaged (in time and space) process. 

This means that z^'(-R) = v"{R) = at some R = Rc- These two equations and the 
dependence v{R) itself, not only determine -Rc and v^. = z/(i?c) but also provide an additional 
relation that involve other parameters of the problem which clearly enter the function i^(i?). 
These are the Mach number M, the heating rate a and the maximum momentum Pmax = 



Pi. For example, given M and a we can easily calculate Pmax 122|. These results show, 
that at least in shocks with high Mach numbers and no significant turbulent heating, the 
nonlinear effects may significantly limit the maximum energy achievable by this acceleration 
mechanism. Before any relevance of this result to, e.g., the lack of observational evidence of 
TeV protons in SNR can be seriously debated, a number of difficult issues must be addressed 
which we discuss briefly in the next section. 



IV. DISCUSSION 

Instead of relying on parameterization, and guessing about the possible values of rele- 
vant parameters it was suggested that, due to the self-regulation, the system should evolve 
precisely into the critical point within given parameter space [^. A comprehensive evalua- 
tion of this suggestion and the study of self-regulation mechanisms poses a serious challenge 
to plasma physics. Note, however, that the critical self-organization approach was proven 
very useful in describing transport phenomena in laboratory plasmas, when the turbulence 
is generated by transport driving gradients, i.e., the pressure or density gradients, e.g., P5| . 
The most important current issue is the construction of adequate links between the inter- 
nal parameters (and so the Alfvenic turbulence), and their dependence upon the external 
parameters. This is an interdisciplinary problem in high energy astrophysics and nonlinear 
plasma physics. It combines hydrodynamics, particle kinetics with a strong emphasis on 
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the theory of dynamical chaos and ergodicity (injection), Alfvenic turbulence, coUisionless 
nonhnear plasma phenomena and, in particular, the theory of coUisionless shocks. 

Currently not all the aspects of this problem are well understood. We believe that we 
understand how the acceleration operates in the test particle (linear) regime, and have an 
analytic description of nonlinearly modified shocks including particle spectra, although under 
a prescribed Bohm-type particle diffusion. Injection can be calculated currently only for a 
relatively strong subshock, while we need it for varying shock parameters of dynamically 
evolving CR shocks. Turbulence dynamics, associated plasma heating and the turbulent 
transport of CRs are understood to even lesser extent. A systematic determination of 
parameters involved should allow the self-consistent calculation of particle spectra. 

Specifically, this requires detailed understanding of several physical phenomena in a shock 
environment, which include: 

1. generation of waves by energetic protons, wave transport and spectral evolution 

2. turbulent transport of protons, in space and momentum, injection 

3. electron injection triggered by proton generated turbulence 

4. heating by Alfvenic and (magneto-) acoustic turbulence, which controls the global 
shock structure 

The resolution of these issues will lead to the further progress in our understanding of particle 
acceleration and associated emission. 
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FIG. 1. Schematic representation of nonlinear ly accelerating shock. Flow profile with a gradual 
deceleration upstream is also shown at the top. 

FIG. 2. Response of shock structure (bifurcation diagram) to the injection of thermal particles 
at the rate v. The strength of the response is characterized by the pre-compression of the flow in 
the CR shock precursor R = ui/uq (see Fig.l). The flow Mach number M = 150; different curves 
correspond to different values of maximum momentum denoted here as piand normalized to mc. 
For each given v and Pmax one or three (arrows) solutions exist. 

FIG. 3. Bifurcation diagram corresponding to the set of response curves shown in Fig. |2|. Since v 
and Pmax are in reality dynamic rather than control parameters, the response curve moves towards 
the bifurcation curve drawn with the heavy line. The resulting state of the system corresponds to 
the critical (inflection) point on this curve, which can be described as a "self-organized critical" 
(SOC) state. 

FIG. 4. Bifurcation diagrams as in Fig. 2 but for different heating rates a (see text), for Mach 
number M = 150 and the cut-off momentum Pmax — 

10^ 
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